Ab initio study of mechanical and thermal properties of GeTe-based and PbSe-based high-entropy chalcogenides

GeTe-based and PbSe-based high-entropy compounds have outstanding thermoelectric (TE) performance and crucial applications in mid and high temperatures. Recently, the optimization of TE performance of high-entropy compounds has been focused on reducing thermal conductivity by strengthening the phonon scattering process to improve TE performance. We report a first-principles investigation on nine GeTe-based high-entropy chalcogenide solid solutions constituted of eight metallic elements (Ag, Pb, Sb, Bi, Cu, Cd, Mn, and Sn) and 13 PbSe-based high-entropy chalcogenide solid solutions: Pb0.99-ySb0.012SnySe1-2xTexSx (x = 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and y = 0) and Pb0.99-ySb0.012SnySe1-2xTexSx (y = 0.05, 0.1, 0.15, 0.2, 0.25 and x = 0.25). We have investigated the mechanical properties focusing on Debye temperature (ΘD), thermal conductivity (κ), Grüneisen parameter (γα), dominant phonon wavelength (λdom), and melting temperature (Tm). We find that the lattice thermal conductivity is significantly reduced when GeTe is alloyed into the following compositions: Ge0.75Sb0.13Pb0.12Te, Ge0.61Ag0.11Sb0.13Pb0.12Bi0.01Te, and Ge0.61Ag0.11Sb0.13Pb0.12Mn0.05Bi0.01Te. This reduction is due to the mass increase and strain fluctuations. The results also show that Ge0.61Ag0.11Sb0.13Pb0.12Bi0.01Te solid solution has the lowest Young’s modulus (30.362 GPa), bulk and shear moduli (18.626 and 12.359 GPa), average sound velocity (1653.128 m/sec), Debye temperature (151.689 K), lattice thermal conductivity (0.574 W.m–1.K–1), dominant phonon wavelength (0.692 Å), and melting temperature (535.91 K). Moreover, Ge0.61Ag0.11Sb0.13Pb0.12Bi0.01Te has the highest Grüneisen parameter with a reduced and temperature-independent lattice thermal conductivity. The positive correlation between ΘD and κ is revealed. Alloying of PbSe-based high-entropy by Sb, Sn, Te, and S atoms at the Se and Pb sites resulted in much higher shear strains resulted in the reduction of phonon velocity, a reduced ΘD, and a lower lattice thermal conductivity.


Mechanical properties calculations
Ab initio calculation method within the density functional theory (DFT) is a very powerful method these days for mechanical properties calculations.In this section, we followed the method of Hongzhi Yao et al. work 1 to calculate the mechanical properties.A scheme has been applied to calculate the elastic constants of the high entropy solid solutions.Starting with the simple Hook's famous law 2,3 that relates the stress components σi with the strain components εj by the relation: Cij is the elastic constants.From knowing the elastic constants, we calculate the other mechanical properties: compliance tensor Sij, Young's modulus (E), Bulk modulus (K), Shear modulus (G), and Poisson's ratio (η).Here, the Voigt -Reuss -Hill (VRH) 4,5 approximation has been used to derive the above mechanical parameters.According to this approximation, the upper and lower bounds for the structural parameters, such as bulk and shear modulus are respectively given by: The dimensionless Kleinman parameter (ζ) is calculated by using the following formula 6 : The machinability index (µM) can be expressed as follows 7 : The formula of Tian et al. 8 was used to calculate the macro Vicker's hardness parameter HV: = 0.92 (   ) 1.137  0.708 (12) Another semi-empirical formula 9 was used to estimate the micro-hardness (H): Lame's constants (λ, μ) are derived from Young's modulus and Poisson's ratio as following 10 : The elastic anisotropy parameter is characterized by the universal anisotropic index A U 11 : Where the subscripts in GV, GR, KV, and VR are the shear and bulk modulus under Voigt and Reuss estimates respectively.
The percentage of anisotropy in compression and shear are given by 12 : The shear anisotropic factor for the {100} shear planes between the [011] and [010] directions (A1) is given by: The shear anisotropic factor for the {010} shear planes between the [101] and [001] directions (A2) is given by: The shear anisotropic factor for the {001} shear planes between the [110] and [010] directions (A3) is given by:

Debye temperature and thermal conductivity calculations
Debye temperature (ΘD) and the average sound velocity (vm) can be calculated as follows [13][14][15][16] : Here, ρ is the theoretical density of the solid solution model, h, kB, and NA are Planck's constant, Boltzmann constant, and Avogadro's number, respectively.M is the molecular weight and n is the number of atoms in the supercell.vs and vl are the transverse (shear) and longitudinal sound velocities respectively.The compressional (longitudinal) waves and shear waves are estimated by using the values of bulk modulus (K) and shear modulus (G) according to Voigt-Reuss-Hill approach based on the following formulas 17 : The acoustic impedance parameter (Z) value of m0 and nine solid solutions can be estimated by using the following formula 18,19 : Clarke's formula for the minimum thermal conductivity (κmin) is given as follows 20,21 : Cahill's formula for the minimum thermal conductivity is given by 22,23 : Slack's formula 24,25 for thermal conductivity (κ) or for lattice thermal conductivity (κL) (in case of the electronic part of thermal conductivity (κe) is negligible) is given by: Where V is the volume of the supercell, A is constant that can be approximated to be 3.1e-6 25 when κ in W.m -1 .K -1 .γα is the acoustic Grüneisen parameter, Ω is the volume per atom, and T is the temperature in Kelvin unit.
The acoustic Grüneisen constant (γα) was calculated by using the following formula 26 : Another formula that can be used to estimate the Grüneisen parameter 27,28 is: Another formula that can be used to estimate the Poisson's ratio 29 Julian 24,27,30 derived the following formula for the constant A in the Slack's formula( formula (S29)) for κ: The mixed model 31 can give an empirical formula for lattice thermal conductivity(κL) as follows: For complex materials, κL mainly has two contributions: the acoustic (κa) and optical (κo) contributions separately, which are given by 31 : ) (37) An empirical formula developed by Fine et al. 32 is used for calculating a rough estimation of melting temperature(Tmelt) of m0 and nine solid solutions(GeTe-based high-entropy chalcogenide models) with a standard error about ±300 K by using the elastic constants: Thermal expansion coefficient (α) can be roughly estimated from the following formula using the value of shear modulus 33 : The maximum value of phonon wavelength (λdom) can be roughly estimated at any temperature by using the following formula 20 : A different empirical formula was used for calculating an estimation of melting temperature(Tmelt) of Pb0.99-ySb0.012SnySe1-2xTexSx(x=0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and y=0) and Pb0.99-ySb0.012SnySe1-2xTexSx(y=0.00, 0.05, 0.10, 0.15, 0.20, 0.25, and x= 0.25) solid solutions with a standard error about ±300 K by using the elastic constants 32,34,35 :

Fig. S4
The three-dimensional contour plots of the Young's modulus for m0 and other nine solid solutions.

Fig. S5
The three-dimensional contour plots of the linear compressibility (1/K) for m0 and other nine solid solutions.
Fig. S6 the three-dimensional contour plots of shear modulus for m0 and other nine solid solutions.

Fig. S7
The three-dimensional contour plots of the Poisson's ratio for m0 and other nine solid solutions.Table S7.Calculated bulk and shear modulus under Voigt and Reuss estimates (KV, KR, GV, GR), the universal anisotropic index (A U ), and the shear anisotropic factors (A1, A2, A3) for ten models in Ge-Te based high-entropy chalcogenides.

Fig. S1 .
Fig. S1.(a) Ball and stick structure of Ge0.61Ag0.11Sb0.13Pb0.12Bi0.01Te(m5) with their solid solution composition and model numbers in the box at right.(b) Ball and stick structure of Pb0.99Sb0.012Se0.5Te0.25S0.25 with their solid solution composition in the box on the right.

Fig. S18 .
Fig. S18.Gruneisen parameter versus the ratio vl/vs for the ten Ge-Te based high-entropy chalcogenide models.

Fig. S19
Fig. S19Thermal conductivities vs Debye temperature for the ten Ge-Te based high-entropy chalcogenide models.

Table S1 .
Number of atoms for each element of the ten models in GeTe-based high-entropy chalcogenides.

Table S2 .
Calculated lattice parameters of the ten models in GeTe based high-entropy chalcogenides.

Table S5 .
The calculated elastic constants Cij (GPa) for ten models in Ge-Te based high-entropy chalcogenides.

Table S14 .
Calculated bulk and shear modulus under Voigt and Reuss estimates (KV, KR, GV, GR), the universal anisotropic index (A U ), and the percentage of anisotropy in compression and shear factors (Acomp,